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ABSTRACT 

This  paper  consists  of  two  parts.  In  the  first,  a  critical 
analysis  of  well-known  procedures  for  the  computation  of  one¬ 
dimensional  shocked  flows  is  made,  in  order  to  show  the  incon¬ 
veniences  of  computing  finite  differences  across  a  discontinuity  and 
to  prove  that  the  use  of  the  enuations  of  motion  in  conservation  form 
does  not  make  the  results  any  more  accurate.  In  the  second,  a 
technique  is  developed  to  treat  one -dimensional  inviscid  problems 
and  it  is  applied  to  the  problem  of  an  accelerating  piston.  Practical 
and  safe  ways  to  predict  the  formation  of  a  shock  and  to  follow  it 
up  in  its  evolution  are  given. 
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Introduction 

The  numerical  treatment  of  one -dimensional  problems  has  been 
the  object  of  a  great  number  of  papers,  mostly  intended  to  provide  theo¬ 
retical  and  experimental  comparisons  between  different  numerical  tech¬ 
niques.  Such  comparisons  have  failed  to  provide  a  positive  contribution  to 
the  state  of  the  art.  They  have  been  unable, for  example,  to  point  out  why 
the  results  of  the  computations  are  generally  affected  by  spurious  oscilla¬ 
tions,  and  why  their  accuracy  is  so  poor,  although  some  qualitative  agree¬ 
ment  with  exact  solutions,  when  available,  can  be  observed. 

In  Ref.  1  a  general  line  of  attack  for  gas  dynamics  problems  has 
been  presented,  the  philosophy  of  which  is  that  computations  can  be  per¬ 
formed  in  any  number  of  space  dimensions  and  time  with  as  great  an 
accuracy  as  desired.  Such  a  goal  can  be  achieved  not  only  without  increas¬ 
ing  the  computational  time  beyond  unacceptable  limits  but  rather  reducing 
it  by  at  least  one  order  of  magnitude. 

A  test  of  the  technique  outlined  in  Ref.  1  is  given  here  for  one- 
dimensional  problems.  A  detailed  analysis  of  some  of  these  problems  is 
necessary  in  order  to  explain  the  physical  subtlety  of  certain  non-linear 
phenomena  and  of  their  interpretation  from  the  view  point  of  numerical 
analysis . 

In  the  first  part  of  the  paper,  the  problem  of  a  steady  shock  sepa¬ 
rating  two  regions  of  uniform  flow  is  examined  in  an  attempt  to  explain  the 
origin  of  certain  well-known  numerical  difficulties.  In  the  second  and  third 
parts,  an  increasingly  sophisticated  technique  is  introduced,  to  show  how 
such  difficulties  can  be,  one  by  one,  eliminated.  A  demonstration  of  how 
the  technique  can  be  applied  to  very  complicated  one-dimensional  inviscid 
problems  will  follow  in  another  paper. 
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1 .  Shocks  and  the  equations  in  divergence  form 

In  the  typical  one-dimensional  problem  analyzed  in  the  literature, 
a  shock  proceeds,  at  a  constant  speed,  into  a  gas  at  rest.  The  region 
behind  the  shock  is,  then,  uniform.  The  flow  pattern  can  be  imagined  to 
be  produced  by  a  piston  suddenly  set  into  motion  at  a  constant  speed 
(Fig.  1).  The  shock  is  not  assumed  to  be  a  discontinuity;  the  numerical 
techniques  intend  to  provide  a  way  of  representing  the  shock  as  a  fast, 
but  continuous,  transition  between  two  constant  states. 

If  the  approach  outlined  in  Ref.  1  were  used  for  this  problem,  the 
flow  field  would  be  divided  into  two  regions  of  continuous  flow  and  the 
shock  confined  to  the  common  boundary.  Then  the  problem  becomes 
trivial;  in  both  regions  separated  by  the  shock  the  flow  is  uniform,  all 
space  derivatives  vanish  identically  and  consequently  the  time  derivatives 
vanish  identically.  There  is  no  doubt  that,  regardless  of  the  technique 


FIG.  I  FLOW  PRODUCED  BY  A  PISTON  MOVING  AT 
A  CONSTANT  SPEED 
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being  used,  the  two  constant  states  are  maintained.  As  far  as  the  shock 
itself  is  concerned,  a  single  application  of  the  Rankine-Hugoniot  conditions 
provides  the  solution  at  any  instant  of  time  since  the  shock  velocity  is  con¬ 
stant.  In  the  first  part  of  this  paper,  however,  we  will  re-examine  the 
simple  flow  problem  outlined  in  Fig.  1  to  discuss  the  basic  philosophy  of 
the  "smeared  shock"  approach  and  to  justify  our  choice  of  the  other  approach 
as  outlined  in  Ref.  1. 

There  is  no  doubt  that  if  one  could  handle  shocks  by  the  same  tech¬ 
nique  used  in  regions  of  continuous  flow,  the  general  logic  of  a  computa¬ 
tional  program  would  be  highly  simplified,  particularly  when  a  flow  is 
expected  to  be  traversed  by  shocks  in  both  directions  and  when  the  problem 
depends  on  more  than  one  space  variable.  However,  anyone  who  has  some 
experience  in  computing  by  the  method  of  characteristics  knows  that,  as 

soon  as  the  characteristics  coalesce  and  a  shock  builds  up,  the  method 
breaks  down.  The  only  way  the  computation  can  be  continued  consists 
of  fitting  the  shock  as  a  discontinuity,  separating  two  regions  of  continu¬ 
ous  flow  where  the  method  of  characteristics  is  applicable.  The  physical 
nature  of  the  shock  is  reflected  in  the  numerical  procedure  by  the  fact 
that  each  shock  point  is  reached  by  three  characteristics,  two  in  the  super¬ 
sonic  region  and  one  in  the  subsonic  region,  plus  the  particle  path  in  the 
supersonic  region  (Fig.  2)  (the  Mach  numbers  being  computed  from  the 
velocities  relative  to  the  shock).  Therefore,  if  one  pretends  to  ignore  the 
shock,  one  is  forced  to  face  a  redundancy  of  data  which  are  inconsistent  as 
long  as  no  jump  is  permitted. 

It  is  rather  surprising  that  so  much  effort  has  been  spent  against 
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X 

FRAME  RELATIVE  TO  THE  SHOCK 


FIG.  2  CHARACTERISTIC  PATTERN  ON  BOTH  SIDES  OF  A  SHOCK 

the  evidence  offered  by  the  method  of  characteristics  and  even  against  a 
basic  principle  of  mathematics,  as  we  will  see  later  on.  There  is  appar¬ 
ently  only  one  possible  explanation,  which  we  will  try  to  present  here,  as 
plainly  as  possible. 

1 . 1  The  smearing  out  of  a  shock  according  to  von  Neumann  and  Richtmyer. 

Let  us  consider  the  problem  in  its  historical  development.  In  1950, 
von  Neumann  and  Richtmyer,  trying  to  solve  the  problem  of  shocked  flows  by 
numerical  techniques,  found  the  amount  of  labor  involved  too  great,  even 
using  automatic  computers  (in  1950,  programming  for  a  machine  like  ENIAC 
was  extremely  cumbersome;  FORTRAN  was  invented  many  years  later). 

Von  Neumann  and  Richtmyer  had  in  mind  multidimensional  shocked  flows. 

A  treatment  similar  to  the  one  outlined  in  Ref.  1  was,  at  that  time,  beyond 
reach,  as  they  stated:  "Shock  calculations  by  direct  application  of  the  Hugoniot 
equations  would  ordinarily  be  prohibitively  difficult,  even  for  rapid,  auto¬ 
matic  computers."  (Ref.  2)  Consequently,  they  suggested  a  device, inspired 


by  the  physical  behavior  of  a  gas  crossing  a  shock.  When  viscosity  and 
heat  conduction  are  taken  into  account,  shocks  are  naturally  smeared 
out  and  can  be  described  as  a  continuous  transition.  The  larger  the  vis¬ 
cous  effects,  the  wider  the  transitional  region.  It  is  known  that,  except 
when  the  gas  is  rarefied,  this  region  is  still  too  thin  to  be  perceivable. 
However,  von  Neumann  and  Richtmyer  thought  that  an  artificial  viscosity 
could  be  forced  into  the  equations  of  motion  in  such  a  way  that  (i)  the 
shock  transition  would  be  spread  out  over,  say,  a  couple  of  mesh  intervals 
(which  meant  to  have  an  unrealistically  large  viscosity  in  that  region),  but 
also  (ii)  managing  to  keep  the  viscous  effects  negligible  in  the  rest  of  the 

computational  mesh.  The  suggested  artificial  viscosity  for  one-dimen- 

2 

sional  problems  depended  on  a  term  containing  u^. 

The  aim  of  the  suggested  device  is  to  compute  inviscid  flows  with 
numerical  methods  apt  to  approximate  continuous  flows,  without  having  to 
make  special  provisions  for  shocks  in  advance,  but  reaching  a  good  approxi¬ 
mation  to  the  inviscid  flow  solution,  except  on  a  two  mesh  interval  bracket¬ 
ing  each  shock. 

1 . 2  Lax's  technique;  artificial  viscosity  and  the  equations  in 
conservation  form. 

In  1954,  Lax  (Ref.  3)  suggested  two  modifications  to  the  above 
approach.  The  first  had  to  do  with  a  way  of  introducing  artificial  viscosity. 
To  see  the  idea  without  complications,  consider  the  equation 

£«  ■  Sx  <» 

Let  the  points  to  be  computed  (  nodal  points)  be  equally  spaced  along  the 

x-axis,  so  that  their  abscissae  are  x  =  nAx  (n  =  1,2,  3, .  .  .  . ).  Let  k 

denote  a  value  of  time  t,  =  kAt  .  Let  f^  denote  the  value  of  f  at  x  =  x 

k  n  n 

2 ) / 2  Ax 


and  t  = 


V 


Approximate  gx  by 


centered  differences  (gn+j  -  g 


’n- 


I 


(i 


but,  instead  of  writing 

k  k 

gn+l  "  gn-l  (2) 

+  - ZTx - 

as  one  would  do  by  a  straight  forward  application  of  Euler's  integration 
rule,  let  us  write 


^k+1  _  ^k 
n  ~  n 


.k  Jc  k  k 

jk+1  _  n+1  n-1  +  gn+l  gn-l  ^ 

2  2  Ax 


n 


(3) 


In  other  words,  the  initial  value  of  f  used  to  approximate  the  time  deriva¬ 
tive  at  the  n-th  point  is  not  the  value  of  f  at  the  n-th  point  but  the  average 
of  the  values  of  f  at  the  two  neighboring  points. 

On  the  other  hand,  if  we  consider  the  equation 

f  =  g  +  vf  (4) 

t  °x  xx  ' 


with 


(5) 


approximate  f^  by  (fj^+j  +  fj^_j  -  2  f^)  /  Ax^,  and  use  an  integration 
scheme  similar  to  (2),  we  obtain  (3)  again.  Terms  containing  second  order 
spatial  derivatives,  such  as  vf^,  appear  in  the  Navier-Stokes  equations, 
where  v  is  proportional  to  the  kinematic  viscosity  of  the  gas  (in  non-dimen¬ 
sional  form,  v  is  the  inverse  of  a  Reynolds  number).  Therefore,  we  may 
consider  (3)  as  a  formula  solving  a  problem  of  motion,  somehow  affected  by 
a  pseudo-viscosity  whose  Reynolds  number  is  of  the  order  of  1/v. 

The  second  innovation  in  Lax's  technique  consisted  of  approximating 
by  a  finite -difference  scheme  not  the  well-known  Euler's  equations 

°t  +  u°x  +  0  ux  =  0 
put  +  0uux  t  px  = 

st+usx=,° 


0 


(6) 
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but  a  system  written  in  what  is  called  a  divergence  form,  or  conservation 
form. 

Such  a  system  can  be  obtained  by  algebraic  manipulation  of  (6)  or 
directly  from  the  equations  expressing  conservation  of  mass,  momentum 
and  energy.  Let 


E 

H 

q 


P  +  DU 


2 


(Y-l)  E  +  ill 


(7) 

(8) 

(9) 

(10) 


The  conservation  equations  mentioned  above  can  be  written  in  the  form: 


p  =  -  m 
^t  x 


mt  =  -  9x 


E  =  -  H 
t  x 


(11) 


and  it  is  easy  to  see  that  (11)  is  equivalent  to  (6).  Formally,  the  system 
consists  now  of  linear  relations  between  first  derivatives.  The  system,  how¬ 
ever,  is  obviously  non  linear  since  the  dependent  variables  are  related  to  one 
another  by  the  non  linear  equations  (7),  (9),  and  (10). 

Note  that  (1)  may  represent  all  three  equations  of  motion  if  f  and  g 
are  considered  as  vectors,  defined  by 


f  = 


D 

ml 

E 


g  =  - 


m 

q 

H 


(12) 


What  happens  if  the  motion  is  steady  is  remarkable  (and,  incidentally, 
the  arguments  and  conclusions  which  follow  are  not  limited  to  one -dimensional 
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problems).  The  left-hand  sides  of  the  equations  vanish  identically.  There¬ 
fore,  the  same  is  true  of  the  right-hand  sides.  And,  of  course,  the  state¬ 
ment  holds  wherever  the  original  equations  (11)  are  applicable,  that  is, 
throughout  all  region  where  the  flow  parameters  are  continuous  and  differ¬ 
entiable.  On  the  other  hand,  consider  the  Rankine-Hugoniot  conditions 
across  a  shock , 

Pi  Vj  =  p2  v2 

P!  vx2  +  Pj  =  P2  v22  +  p2  (13) 

,  .  1  2  ,  ,  1  2 
hl+7vl  =  h2  +  Iv2 

where  v  is  the  velocity  relative  to  the  shock,  h  is  the  enthalpy: 

h  =  -ly  £  (14) 

Y-*  P 

and  the  subscripts  1  and  2  refer  to  conditions  on  either  side  of  the  shock. 

In  a  steacy  state,  the  velocity  v  coincides  with  u.  By  using  the  definitions 
(7),  (10)  and  (9),  the  Rankine-Hugoniot  conditions  can  be  written  as 

ml  =  m2 

qL  =  q2  (15) 

Hi  =  H2 

In  other  words,  the  quantities  m,  q,  and  H  (in  physical  terms ,  mass  flow, 
momentum  and  total  enthalpy)  are  continuous  across  a  shock.  And  since 
m  =  0,  q  =  0,  H  =  0  everywhere  else,  m,  q,  and  H  are  constant 
throughout  the  flow.  It  is  obvious,  then,  that  any  finite-difference  approxi¬ 
mation  to  the  x-derivatives  of  m,  q,  and  H  vanishes  identically,  even  if  the 
nodes  lie  at  opposite  sides  of  a  shock.  By  using  the  equations  in  the  form  (11) 
their  finite -difference  counterparts  can  be  used  all  throughout  a  shocked  flow 
even  when  the  individual  components  of  mass,  momentum  and  total  enthalpy, 
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that  is,  density,  velocity  and  pressure,  are  not  continuous  (and  thus  not 
differentiable)  across  a  shock.  Let  us  hasten  to  say  that  a  statement  of 
this  kind,  although  obviously  true,  is  a  deceiving  one.  What  it  actually 
means  is  that  if  we  compute  a  steady  flow  by  finite -difference  technique 
as  if  it  were  unsteady,  and  we  use  as  initial  conditions  the  correct  steady 
flow,  and  we  use  the  equations  in  conversation  form,  we  find  that  all  time 
derivatives  are  zero,  and  therefore  we  compute  nothing  at  all.  The  result 
is  self-consistent  but  is  trivial. 

It  is  surely  a  better  result  than  what  we  would  have  got,  had  we 
used  a  finite  difference  counterpart  to  (6)  in  the  presence  of  a  shock  since 
finite-difference  approximations  to  the  space  derivatives  across  a  steady 
shock  would  generally  provide  finite  values  of  the  time  derivatives,  and 
therefore  local  changes  in  the  physical  parameters,  a  result  which  is 
inconsistent  with  the  steady  state  assumption. 

The  question  is,  Does  then  a  set  of  finite -difference  equations 
obtained  from  (11)  give  better  results  than  a  set  obtained  from  (6)  when 
the  flow  is  not  steady?  Lax  himself  did  not  attempt  to  prove  anything  of 
the  sort;  he  only  advanced  some  hopeful  conjectures.  Then  he  made  some 
numerical  work  to  try  out  his  innovations.  Fig.  3  is  a  plot  of  the  pres¬ 
sure  distribution  as  described  in  Table  V  of  Lax's  paper,  (Ref.  3).  The 
theoretical  distribution  is  a  constant  line,  at  p  =  50,  followed  by  a  con¬ 
stant  line  at  p  =  0.  The  transition  should  occur  approximately  where  indi¬ 
cated.  The  figure  shows  a  typical  "qualitatively  good"  result.  Let  us 
analyze  it,  however,  keeping  in  mind  that  two  questions  must  be  answered,  viz. 

(i)  Is  the  transition  replacing  the  shock  sufficiently  sharp?,  and  (ii) 

Does  the  conservation  form  of  the  equations  solve  the  problem  of  computing 
across  a  shock? 

Two  features  of  the  figure  should  be  considered  carefully.  The 
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first  is  the  number  of  mesh  intervals  along  which  the  shock  is  spread. 

The  second  is  a  mild  oscillation  appearing  in  the  high  pressure  side  of  the 
"shock"  (to  make  it  more  evident,  part  of  the  plot  has  been  redrawn  at  a 
larger  scale  in  the  upper  part  of  the  figure).  These  two  features  are  typ¬ 
ical  of  numerical  effects  which  are  unrelatable  to  the  physical  nature  of 
the  flow.  We  are  going  to  discuss  them  in  detail  in  the  first  part  of  this 
paper,  showing  that  the  oscillation  is  the  typical  symptom  of  numerical 
trouble  in  attempting  to  deal  with  a  discontinuity  by  "continuous"  tech¬ 
niques,  whereas  the  spreading  of  the  transition  over  too  many  mesh 
points  is  the  telltale  of  a  gigantic  artificial  viscosity. 

Unfortunately,  in  Lax's  procedure  one  effect  cannot  be  studied  with- 

$ 

out  having  the  other  as  well,  and  conclusions  cannot  be  easily  drawn.  How¬ 
ever,  one  can  grasp  the  idea  by  comparing  Figs.  3  and  4  where  some  results 
obtained  by  Emery  (Ref.  4)  with  the  Lax  scheme  are  replotted.  The  curve  at 
the  left  and  the  one  at  the  right  of  Fig.  4  were  obtained  by  using  At/ Ax  =  .9  and 
At/Ax  =  .7,  respectively.  The  curve  of  Fig.  3  corresponds  to  At/ Ax  =  .25. 

In  a  real  shock,  for  a  given  Mach  number  (not  much  higher  than  4),  the 
shock  thickness  5  is  roughly  proportional  to  the  inverse  of  the  Reynolds 

number.  In  Lax's  scheme, (5)  shows  that  the  Reynolds  number  is  propor- 
2  2 

tional  to  At/Ax  .  Therefore,  the  product  5(At/Ax  )  should  be  practically 
the  same  for  all  three  cases  if  the  artificial  viscosity  acts  as  a  real  viscos¬ 
ity.  Now, 6 /Ax  is  the  number  of  intervals  over  which  the  "shock"  is  spread. 


This  is  due  to  the  fact  that  the  Euler's  integration  rule  (2)  is  uncondi¬ 
tionally  unstable  and  (3),  as  we  have  seen,  contains  artificial  viscosity. 
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From  the  figures  we  find 


6  /  Ax 

At/ Ax 

6  (At/ Ax2) 

12 

.25 

3 

5 

.7 

3.5 

4 

.9 

3.6 

From  the  figures  of  the  last  column  we  can  conclude  that  the 
effects  of  artificial  viscosity  on  shock  thickness  are  very  similar  to  the 
effects  of  real  viscosity.  Unfortunately,  we  can  also  see  that  the  artificial 
Reynolds  number  per  unit  length  is  of  the  order  of  1/Ax.  In  practice  one 
deals  with  intervals  which  cannot  be  smaller  than,  say,  1/50  of  a  typical 
reference  length,  lest  the  computational  time  increases  beyond  acceptable 
limits.  Therefore  the  artificial  Reynolds  number  is,  at  most,  of  the  order 
of  50. 

Figs.  3  and  4  show  that  the  artificial  viscosity  plays  a  crucial  role 
in  damping  out  the  oscillations  and  that  an  effective  damping  is  obtained 
only  if  the  artificial  Reynolds  number  is  extremely  low,  so  that  the  transi¬ 
tion  replacing  the  shock  cannot  be  sufficiently  sharp. 

1.3.  The  Lax-Wendroff  Technique. 

Reducing  the  artificial  viscosity  in  Lax's  scheme  entails  increasing 
the  amplitude  of  the  oscillations  in  the  high-pressure  side  of  the  shock. 
Therefore  we  cannot  evaluate  whether  the  conservation  form  of  the  equation 
of  motion  is  of  any  help  in  computing  a  flow  which  contains  a  sharp  disconti¬ 
nuity.  For  such  an  analysis,  one  needs  a  scheme  which  can  work  in  the 
absence  of  artificial  viscosity.  The  scheme  suggested  by  Lax  and  Wendroff 
in  I960  (Ref.  5)  seems  to  be  aimed  at  increasing  the  accuracy,  and  eliminating 
artificial  viscosity,  while  maintaining  the  formulation  of  the  equations  in  con¬ 
servation  form,  obviously  in  the  hope  that  shocked  flows  could  be  treated 
without  shock  fitting. 


The  Lax-Wendroff  scheme,  in  principle,  depends  on  approximating 
the  unknowns  at  a  time  t  +  At  by  a  Taylor  expansion  truncated  at  the  second 
order  terms: 

f(t+  At)  =  f(t)  +  ftAt  +  ^fttAt2  (16) 

expressing  the  time  derivatives  through  space  derivatives,  making  use  of 
the  equations  of  motion,  and  replacing  the  space  derivatives  by  centered 
differences.  Lax  and  Wendroff  have  applied  the  scheme  to  the  equations 
of  motion  formulated  in  conservation  form.  In  the  one -dimensional  case, 
(11)  is  the  starting  point.  Then  (16)  is  used.  To  compute  ftt,  (1)  is 
differentiated  with  respect  to  t: 

£tt  =  gtx 
but 


gt  =  J  ft  =  J  gx 

(where  J  is  the  matrix  whose  determinant  is  the  Jacobian  of  g  with  respect 
to  f),  so  that 


£tt  =  <J  gx>x 


(17) 


For  no  apparent  reason,  (J  gx)x  is  not  resolved  into 


<J  «x>x  =  Jxgx  +  J  gxx 


but  is  computed  as 


<J  8x>x  = 


I  0k  k  _k  .  _k 

.  |Jn+l+Jn  6n+1-gn  J-  +J. 

2 


n 


n-l  *n  -  «n-l 


Ax 


(18) 


(19) 


Ax  2  Ax 

which  implies  an  averaging  of  J  and,  thus,  a  higher  truncation  error  than 
if  (18 )  were  used. 

We  would  like  to  mention  in  passing  that  the  scheme  suggested  by  (16) 


does  not  require  the  equations  of  motion  to  be  cast  in  conservation  form. 
From  a  practical  point  of  view,  it  is  more  convenient  to  use  Euler's 
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equations.  Recasting  the  equations  of  motion  in  conservation  form  is  a 
cumbersome  process,  particularly  when  more  than  one  space  dimension  is 
involved.  Auxiliary  mappings  of  the  physical  space  onto  a  computational 
space  make  the  procedure  even  more  complicated,  if  not  impossible.  The 
Jacobian  matrices  in  multi-dimensional  problems  not  only  are  complicated 
but  require  a  great  amount  of  numerical  computations  which  is  avoided  if 
Euler's  equations  are  used.  In  conclusion,  a  program  based  on  the  original 
Lax-Wendroff  suggestion  is  more  complicated  (and  consequently  more 
exposed  lo  accidental  mistakes),  much  more  time  consuming  on  the  machine 
and  less  accurate  than  a  program  applying  (16)  to  Euler's  equations. 

Of  course,  these  disadvantages  would  be  largely  balanced  by  the 
possibility  of  treating  shocks  as  sharp  transitions .  However,  the  equations 
in  conservation  form  do  not  keep  up  to  the  hopes  suggested  by  Lax's  con¬ 
jectures.  Fig.  5  shows  a  typical  result  obtained  by  Lax  and  Wendroff.  Here 
the  oscillations  in  the  high-pressure  side  of  the  shock  are  much  stronger 
than  in  Fig.  1.  The  figure  is  borrowed  from  page  333  of  Richtmyer  and 


FIG.  5  VELOCITY  DISTRIBUTION  ACROSS  A 
SHOCK,  AS  COMPUTED  BY  LAX  AND 
WENDROFF 


lb 


Morton's  book  (Ref.  6),  where  no  explanation  is  given  for  the  appearance  of 
the  oscillations.  It  is  only  stated  that  the  oscillations  can  be  reduced  by 
adding  again  an  artificial  viscosity  to  the  equation  of  motion. 

1 . 4  Explaining  the  formation  of  wiggles. 

From  the  examolcs  above  it  appears  that  wiggles  tend  to  form  in  the 
high-pressure  side  of  a  shock,  regardless  of  the  computational  technique 
being  used,  unless  artificial  damping  devices  are  used.  Here  we  attempt  to 
explain  the  formation  of  such  wiggles.  We  are  interested,  of  course,  in  a 
numerical  scheme  whose  artificial  viscosity,  if  any,  is  such  that  the  shock 
is  not  spread  over  more  than  one  mesh  interval.  We  consider  a  shock 
moving  from  left  to  right  at  a  constant  speed,  W,  and  separating  two  regions 

of  uniform  flow.  Consider  three  points,  labeled  n-l,n,  and  n+1  at  the 
abscissae  x-Ax,  x,  and  x+Ax  (Fig.  6)  and  suppose  that,  when  a  certain 
computation  is  started,  the  shock  is  somewhere  between  x-Ax  and  x. 

Even  if  we  are  willing  to  accept  a  continuous  transition  instead  of  the 
jump  at  the  shock,  we  must  acknowledge  that  the  numerical  scheme  is 


P 


FIG.  6  DISCONTINUITIES  ACROSS 
A  SHOCK 
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unable  to  distinguish  between  such  a  transition  and  a  jump, since  both  take 
place  between  two  adjacent  nodes  and  the  values  of  the  physical  parameters 
are  defined  only  at  nodes.  Under  these  conditions,  we  can  show  how 
oscillations  build  up  long  before  a  steady  state,  when  physically  existing, 
can  be  reached. 

In  the  above  problem,  the  three  quantities  m,  q,  and  H  are 
higher  at  the  left  of  the  shock  than  at  the  right  .  The  same  may  be  said 
of  P,  p,  u,  and  E.  Let  F  be  any  of  these  quantities  and  let  Fig.  7a 
represent  what  the  F-distribution  should  actually  be  at  a  certain  time  t. 

In  the  same  figure, a  set  of  initial  values  for  a  numerical  computation  is 
also  shown,  at  a  discreet  number  of  nodal  points.  At  all  nodal  points 
except  A  and  B  the  results  of  a  numerical  computation  are  correct.  We 
are  going  to  analyze  the  results  of  one  computational  step  at  points  A  and 
B  by  assuming  that 


Such  a  condition  can  always  be  satisfied  if  At  is  chosen  sufficiently  small. 
Note  then  that  the  numerical  approximations  to  the  derivatives  at  A  and  B, 
as  shown  in  Fig.  7b  by  the  dotted  lines,  are  negative;  according  to  (11),  all 
quantities  F  increase  at  A  and  B.  At  t  +  At,  the  actual  F-distribution 
(where  the  shock  moved  somewhat  ahead)  and  the  computed  F-distribution 
are  as  in  Fig.  7c.  Now  points  C  and  D  start  being  affected  by  errors  and 
we  analyze  the  results  of  a  second  computational  step  with  (20)  still  valid. 
The  approximate  derivatives  are  negative  at  A,  B  and  D,  and  positive  at  C; 
consequently  F  increases  at  A,  B  and  D  and  decreases  at  C  (Fig.  7d. ) 
Without  pursuing  the  argument  any  further,  one  can  already  see  the 


*  They  are  the  same  only  if  the  shock  is  steady. 
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FIG.  7  FORMATION  OF  WIGGLES 


formation  of  a  wiggle  on  the  high  pressure  side  of  the  shock.  This  is  a 
common,  typical  situation,  due  to  a  bold,  but  hardly  justifiable,  attempt 
to  compute  derivatives  across  a  discontinuity;  and  it  is  a  first-order 
effect.  The  introduction  of  higher  order  terms,  as  in  the  Lax-Wendroff 
technique,  is  irrelevant  as  far  as  this  effect  is  concerned. 

However,  as  the  computation  proceeds,  the  second  order  terms  in 
(16)  become  more  and  more  important  and  they  are  finally  responsible  for 
the  partially  oscillatory  and  unsymmetrical  distributions  typified  by  Fig.  5. 

A  simple  example  is  offered  by  the  density  distribution,  which  is  governed 
by  the  first  of  (1  i).  The  second  order  term  in  the  finite-difference  formula¬ 
tion  has  the  sign  of  m  and  it  grows  with  m  .  As  the  amplitude  of  the 

XX  XX 
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1  2 

oscillations  increases,  the  contribution  of  At  ends  by  becoming  of 

the  same  order  as  m^At.  At  point  A  of  Fig.  7d  the  first  order  term  tends 
to  increase  n  but  the  second  order  term,  which  is  negative,  tends  to 
decrease  it.  The  same  can  be  said  of  point  C,  with  the  signs  reversed. 
Consequently,  the  density  distribution  tends  to  freeze  in  an  oscillatory 
pattern. 

In  conclusion,  the  fact  that  a  certain  numerical  technique  has  been 
proven  to  be  a  valid  approximation  of  a  continuous  differentiable  flow  does 
not  make  its  use  legitimate  across  a  discontinuity.  All  mathematical  proofs 
founded  on  differentiability  are  invalid  and  the  numerical  consequences  of 
doing  something  mathematically  wrong  appear  as  wiggles  on  the  high  pres¬ 
sure  side  of  the  discontinuity.  An  accurate  steady  state  solution  cannot  be 
reached  asymptotically  since,  once  the  wiggles  are  formed,  the  numerical 
evolution  has  very  little  in  common  with  the  physical  evolution. 

In  addition,  one  can  easily  see  experimentally  that  a  steady  state  with 
discontinuities  is  numerically  unstable,  if  treated  with  the  Lax-Wendroff 


FIG.  8 

VELOCITY  DISTRIBUTION 

ACROSS  A  STEADY 
SHOCK,  COMPUTED  BY 
THE  LAX-WENDROFF 
SCHEME  WITH  NO 
ARTIFICIAL  VISCOSITY 
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scheme  ,  even  if  the  equations  are  in  conservation  form,  when  the  discon¬ 
tinuities  lie  inside  the  computational  region.  A  minor  perturbation  in  the 
assigned  data  generates  non-zero  derivatives  and  triggers  an  unsteady,  non¬ 
physical  evolution.  See  Fig.  8  where  u  is  plotted  after  computing  20  steps 
of  a  steady  flow  (corresponding  to  a  Mach  2  shock)  by  the  original  Lax- 

Wendroff  technique,  with  At/ Ax  =  .2.  The  initial  values  were  all  exact 

-4  s'{ 

except  that  u,  at  one  point,  had  an  error  of  the  order  of  10 

1 .  5  Inconveniences  of  artificial  viscosity. 

The  argument  developed  in  the  preceding  section  is  based  on  the 
assumption  that  there  is  no  artificial  viscosity,  but  it  holds  as  long  as  the 
artificial  viscosity  is  too  small.  It  is  clear  that  the  shock  transition  cannot 
be  confined  to  a  single  mesh  interval.  It  is  shown  in  Ref.  7  that  the  same 
conclusion  is  reached  for  a  viscous  flow.  Qualitatively,  our  argument  agrees 
with  Crocco's  analysis  of  a  shock  in  a  viscous  flow  (Ref.  11).  In  Crocco's 
paper  the  equations  are  somewhat  simplified  and  a  different  numerical 
scheme  is  used. 

By  so  doing,  it  is  found  that  certain  distributions  of  wiggles  are 
the  only  possible  numerical  solution  of  the  shock  problem  if  the  mesh  size 
is  too  small  with  respect  to  the  shock  thickness.  A  typical  condition  found 
by  Crocco  is,  for  example, 

Ax  <  6/2 

where  6  is  the  shock  thickness-  defined  by  the  ratio  of  the  velocity  jump  to 
the  maximum  slope  of  the  velocity  distribution.  If  Lax-Wendroff 1  s  scheme 

*  The  program  used  to  compute  this  case  is  reported  in  Appendix  1 .  For  the 
reader's  peace  of  mind  it  should  be  mentioned  that  no  instabilities  appear  if 
the  initial  values  are  all  exact. 
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is  used,  the  wiggle  distribution  about  the  shock  is  asymmetrical.  The 
amplitude  of  the  oscillations  is  much  higher  on  the  high-pressure  side  than 
on  the  low-pressure  side. 

The  above  considerations  show  that  the  perspectives  for  a  practical 
usage  of  artificial  viscosity  are  rather  gloomy.  Of  course,  artificial  viscosity 
makes  the  numerical  procedure  once  more  mathematically  legitimate. 

However,  one  faces  the  unpleasant  choice  between 

(i)  a  large  viscosity  (to  maintain  a  coarse  mesh  and  a  short  computa  ¬ 
tional  time);  but,  then,  the  problem  solved  numerically  is  not  the  original 
inviscid  problem. 

(ii)  a  fine  mesh  (to  reduce  the  shock  thickness  and  the  dissipative 
effects  in  general  to  acceptable  values);  but  then,  the  price  to  be  paid 

in  terms  of  computational  time  becomes  too  high. 

In  the  simple  problem  examined  in  the  preceding  sections,  the  dam¬ 
aging  effects  of  artificial  viscosity  are  not  as  strong  as  in  more  complicated 
problems.  In  fact,  in  most  of  the  computational  region  the  flow  is  uniform; 
all  derivatives  are  zero,  and  the  artificially  viscous  terms  are  zero.  In 
multi-dimensicnal  problems  where  non-uniform  steady  states  can  be  reached, 
the  steady  state  affected  by  artificial  viscosity  is  generally  different  from  the 
physical  one  for  inviscid  flow.  For  example,  if  Lax's  technique  is  used,  the 
values  of  all  physical  parameters  become  so  distributed  that  their  time  incre¬ 
ments  exactly  balance  the  errors  due  to  averaging  (Ref.  8). 
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2.  The  problem  of  the  accelerating  piston. 


We  are  going  now  to  examine  a  problem  which  is  more  interesting 
than  the  shock  separating  two  uniform  regions  and  which  leads  to  a  number 
of  practical  applications.  It  is  the  problem  of  what  flow  is  produced  in  a 
gas  at  rest  by  a  piston  moving  at  variable  speed.  Let  x  =  b(t)  be  the  tra¬ 
jectory  of  the  piston.  Let  us  assume  that  the  piston  was  kept  at  rest  when 
t  <  0  and  its  motion  starts  at  t  =  0.  Let  us  also  assume  that  the  initial 
speed  of  the  piston,  b(0),  is  zero.  For  t  >  0  we  will  assume  that  b(t)  is 
analytic  and  regular. 

Let  the  pressure  and  density  of  the  gas  at  rest  be  the  reference 


values,  Pre£  and  Pref»  and  let  a  reference  length,  x^^,  be  chosen  arbi¬ 
trarily.  A  reference  velocity  is  then  defined  as  v'Pj.gf/  pref  and  a  re^er" 
ence  time  as  xref/ °ref  •  addition,  let 


7  =  £  fjL 2*  ,  P  =  tn-£— ,  S  =  y  In  7  -  (Y  -  U  P  (21) 

0  Pref  Pref 


In  non-dimensional  form,  the  equation  of  motion  can  be  written  as  follows  : 

P.  +  uP  +  yu  =  0 

t  X  X 

u.  +  uu  +  ?  P  +  0  (22) 

t  x  x  '  ' 


St  +  uS  =  0 
t  x 


The  characteristics  are  defined  by 

\  -  -  u  +  a  (23) 

where  a  is  the  local  speed  of  sound,  a  =  -Jy  .  We  will  say  that  the 
characteristics  belong  to  the  first  or  to  the  second  family  if  the  upper  or  lower 


#  The  theoretical  background  supporting  the  information  contained  in  equations 
(28)  to  (31)  is  assumed  as  known  ar.d  can  be  found  in  Ref.  9. 


sign  in  (23)  is  taken,  respectively.  The  compatibility  equations  along  the 
characteristics  are 


I  a  dP  +  y  du  =  0 

II  a  dP  -  y  du  =  0 


(24) 


If  the  entropy  is  constant  throughout,  the  compatibility  equations  can  be 
integrated  and  recz.st  in  the  form: 


I  2  ,  , 

— a  +  u  =  2r 

Y  (25) 

II  a  -  u  =  2s 

where  r  and  s  are  constant  along  the  characteristics.  The  first  character¬ 
istic  of  the  I  family  issuing  from  the  point  (x  =  0,t  =0)  is  the  line 

x  =  fyt  (26) 

This  line  divides  the  piston-driven  flow  from  the  gas  at  rest.  As  long  as 
characteristics  of  the  I  family  do  not  coalesce,  the  region  above  such  a 
line  is  a  simple,  isentropically  compressive  wave,  corresponding  to  a  con¬ 
stant  value  of  s  .  The  value  of  u  at  a  point  (x,t)  in  the  wave  can  be  exactly 
computed  as  follows . 

Let  t  be  the  value  of  t  at  which  the  I- characteristic  passing  through 
the  point  issues  from  the  piston  (Fig.  9).  Along  such  line,  u,  a  and  X.  are 
constant,  so  that  the  line  is  straight  and  defined  by 

x  =  b(T)  +  [b(T)  +  aj(T)J(t  -  t)  (27) 

where  a^(T)  is  the  speed  of  sound  on  the  piston  at  time  T.  Since  s  has  the 
same  value  inside  the  simple  wave  and  in  the  gas  at  rest, 
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from  (25.11)  we  obtain 

aj(T)  -  b(T)  =  -^f  (29) 

that  is , 

aj  (T )  =  77  +  ^  b(T)  (30) 

By  replacing  this  value  of  a^(T)  into  (27) 

X  =  b(T)  +  £/y  +  b(T)J  (t  -  T)  (31) 

is  the  equation  of  a  I-characteristic  in  the  simple  wave. 

Suppose  now  that  one  wants  to  know  the  u- distribution  along  the 
x-axis  at  a  given  time  t ^ .  Equation  (31)  should  be  solved  for  t,  after 
replacing  t  by  tj .  Then, 

u(x,  tj)  =  b(T)  (32) 

The  derivative  of  u  with  respect  to  x  is 

ux(x,  t)  =  b(T)  |1  (33) 

where  is  computed  from  (31)  at  constant  t: 


FIG.  9  FLOW  PRODUCED  BY  AN 
ACCELERATING  PISTON 
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dT  _  _ 1 _ 

^  ^b(T)(t-T)-77-X~L  b(77 


(34) 


so  that 


u 


2  b(T) 


(Y+l)  b(T)  (t-T)  -  277-  (Y-l)b(T) 
Similarly,  the  derivative  of  u  with  respect  to  t  is  computed: 


(35) 


ut*“ 


b(T)[277+(Y+l)b(T)] 


(36) 


(V+l)b(T)(t-T)  -  277-  (Y-l)b(T) 

As  long  as  the  denominator  in  (35)  and  (36)  does  not  vanish,  the  derivatives 
of  u  are  continuous  and  differentiable,  except  along  the  first  characteristic, 
defined  by  (26).  if  the  initial  acceleration  of  the  piston  is  not  zero. 

When  the  characteristics  coalesce,  they  define  an  envelope,  which 
can  be  found  by  differentiating  (31)  with  respect  to  T: 

b(T)  (t-T)  -  77  -  b(T)  =  0  (37) 

and  eliminating  T  between  (31 )  and  (37 ),  or  else,  considering  (31 )  and  (37) 
as  the  parametric  equations  of  the  envelope: 


X  =  b(T)  4p/7  +  Iilb(T)|2VvMV-l)b(T) 
L  J  (Y+l)b(T) 


(38) 


t  =  t  +  2  77  +  (Y-l)b(T) 


(Y+Ub(T) 

As  for  the  point  where  the  characteristics  coalesce,  we  must  distinguish 
between  two  cases. 

1)  The  initial  acceleration  of  the  piston,  b(0),  is  different  from 
zero.  In  this  case,  we  can  prove  that  the  curve  defined  by  (  38)  is  tangent 
to  the  first  I-characteristic,  defined  by  (26).  Indeed,  if  t=0  (and,  then, 
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b(T)=0,  b(T)  =  0),  it  is  eas  ily  seen  that  x  and  t  as  defined  by  (38)  satisfy 
(26).  In  addition,  dx/dt,  computed  from  (38),  turns  out  to  be  equal  to  Jy  . 

The  location  of  the  first  characteristic  and  of  the  envelope  is  shown  in  Fig.  10. 
Point  A,  defined  by 


=  2  vV 

(y+l)b(0) 


(39) 


is  their  common  point  and  is  also  the  point  where  the  characteristics  coalesce, 
that  is,  the  origin  of  a  shock  wave. 

It  is  worth  studying  how  the  distribution  of  a  typical  parameter,  say  u, 
steepens  up  in  the  vicinity  of  the  first  characteristic  as  time  increases,  until 
its  slope  becomes  imfinitely  large  at  A.  From  (35)  v'ith  T  =  0,  it  follows  that 


u 


x 


2  b  (0) 

(Y+l)b(0)t  -  2/7 


(40) 


FIG.  10  ENVELOPE  OF  CHARACTERISTICS 
FOR  THE  CASE  OF  FINITE  INITIAL 
ACCELERATION 
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Fig.  11  shows  u  as  a  function  of  t  along  the  first  characteristic,  in  the 

II 

case  where  b(0)  =  2,  y  =  1.4.  Note  that  the  motion  of  the  piston  does  not 
affect  this  steepening  process.  In  fact,  the  right-hand  side  of  (40)  contains 

II 

only  the  parameter  b(0).  Therefore,  regardless  of  what  the  phenomenon 
is  in  the  rest  of  the  simple  wave,  the  value  of  u^  always  becomes  infinite 
at  point  A.  Fig.  12  shows  a  number  of  cross-sectional  plots  of  u(x)  at 
different  times  between  t  =  0  and  t  =  t^  for  the  case  where  the  piston  path 
is  defined  by 

b  =  t2  (41) 


FIG.  II  PLOT  OF  Ux  AS  A  FUNCTION  OF 

TIME  ALONG  THE  FIRST  CHARACTERISTIC 
IN  THE  FLOW  PRODUCED  BY  A  PISTON 
MOVING  WITH  A  FINITE  INITIAL 
ACCELERATION 


*  Drawn  by  a  Stromberg-Carlson  4020  plotter. 
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FIG.  13  PROVING  THE  EXISTENCE  OF  A  CUSP  IN 
THE  ENVELOPE  OF  CHARACTERITICS  WHEN 
THE  PISTON  HAS  NO  INITIAL  ACCELERATION 


2)  The  initial  acceleration  of  the  piston  is  zero.  In  this  case  u^  on 
the  first  I-characteristic,  as  given  by  (40), is  identically  zero  and  no  steep- 
ening  of  the  u(x)  curves  occurs.  This  result  suggests  that  no  coalescence 
of  characteristics  occurs  along  the  first  characteristic.  The  envelope  no 
longer  has  the  form  shown  in  Fig.  10. 

To  find  what  is  the  shape  of  the  envelope,  consider  the  second  of 
(38)  which  gives  t  on  the  envelope  as  a  function  of  T,  At  t  =  0,  t  is  infinite. 

^  *  II 

For  all  T  >  0,  since  b(T)  >  0,  b(T)  >  0  for  an  accelerating  i  ston,  t  is 

1 1 

positive.  If  b(T)  does  not  vanish  for  T  >  0,  then  t  tends  to  »  again  when  T 

1 1 

increases  indefinitely.  If  b(T)  vanishes  at  some  t  =  Tj,  there  t  becomes 
infinite.  In  any  case,  t(T)  on  the  envelope  must  have  a  minimum  (Fig.  13) 
at  a  certain  value  defined  by 


2,/y  +  (Y-l)b(TB) 
b2(TB) 


0 


(42) 


Consequently,  for  increasing  t,  t  first  decreases  and  then  increases,  and 
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FIG.  14  ENVELOPE  OF  CHARACTERISTICS  FOR  THE 
FLOW  PRODUCED  BY  A  PISTON  MOVING 
WITH  THE  LAW  xp*t3 
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the  envelope  is  shaped  as  in  Fig.  14,  with  a  cusp  at  B.  The  location  of  the 

cusp  now  depends  on  the  nature  of  the  piston  path.  Note  that  u^  becomes 

infinite  at  B  because  u^  is  infinite  along  the  entire  envelope,  as  one  can 

see  by  combining  (35)  and  the  second  of  (38).  Therefore,  the  nature  of  the 

u- distribution  at  the  origin  of  a  shock  is  the  same  whether  the  shock  occurs 

at  A  or  at  B.  Fig.  15  shows  a  number  of  cross-sectional  plots  of  u(x)  at 

different  times  between  t  =  0  and  t  =  t0  for  the  case  where  the  piston  path 

Jd 

is  defined  by 

b  =  t3  (43) 


In  the  initial  phase  (prior  to  the  coalescence  of  characteristics),  the 
flow  is  continuous,  although  the  fir-t  derivatives  of  the  flow  variables  may 
be  discontinuous  along  the  I- characteristic  issuing  from  the  origin.  Con¬ 
sequently,  we  may  decide  to  use  a  single  computational  region,  extending 
from  the  piston  (at  the  left)  to  an  arbitrary  boundary  located  somewhere  in 
the  region  of  fluid  at  rest,  where  the  flow  parameters  are  well  known.  For 
simplicity,  assume  a  right  boundary  moving  parallel  to  the  piston  (Fig.  16), 
thus  defined  by 

x  =  c(t)  =  b(t)  +  xq  (44) 

The  computational  region  is  shown  in  Fig.  16;  it  consists  of  a 
segment  of  the  X-axis  and  the  obvious  way  to  relate  the  physical  space  to 
the  computational  space  is  to  let 


X  =  x  -  b(t) 
T  =  t 


(45) 


Since  for  any  function  f, 

f 

x 


(4b) 
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in  the  computational  space  the  equations  of  motion  (22)  are,  in  matrix 


form: 


+  Afx 


0 


with 


■p' 

O 

< 

o 

f  = 

u 

,  A  = 

?  C,  0 

.s_ 

o 

o 

o_ 

(47) 


(48) 


All  interior  points  are  computed  as  follows.  Equation  (47)  is  differentiated 


with  respect  to  X  and  T  successively: 


PISTON 


I 

— 

BOUNDARY 


*o 


X 


FIG.  16  FIRST  DEFINITION  OF  A  COMPUTATIONAL  REGION 
FOR  THE  PROBLEM  OF  THE  ACCELERATING 
PISTON 
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fTX  +  AX  fX  +  A  fXX  =  °  (49) 

fTT  +  At  fx  +  A  fTX  =  0  (50) 

The  first  and  second  space  derivatives  appearing  in  (47  ),  (49)  and  (50)  are 
approximated  by  centered  differences. 


XX 


fk  -  £ 
n+1  n-1 

2  AX 

tl  +  £l  '  2fn 
A3? 


(51) 


(52) 


To  evaluate  and  A,p,  note  that 


(53) 


consequently, 


CX  =  UX*  *x  =  PX  +  7sx>*  ^  = 


cx  0  0 

*x  cx  0 
0  0  c. 


(54) 


These  values  are  used  to  evaluate  f^^.  Similarly, 


0 

0 

C 


(55) 


Tj 


and,  since  u^,,  P^,  and  S^,  are  furnished  by  (47),  we  have  all  the  necessary 
elements  to  compute  f^.  Then  (16)  is  used  to  proceed  one  step  ahead  in 
time.  As  a  matter  of  fact,  the  entropy  is  identically  zero  in  this  problem 
(which  we  consider  finished  once  the  characteristics  coalesce),  the  third 
equation  could  be  dropped  and  the  terms  S,  S^,  S,p  eliminated;  but  we  pre- 
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FIG.  17  TO  DETERMINE  THE  RATIO  BETWEEN 
At  AND  Ax 

fer  to  keep  them  to  show  a  general  pattern  to  be  followed  for  non-isentropic 
flows . 

The  boundary  values  in  the  present  case  are  extremely  easy  to 
compute.  At  the  piston,  u  =  b,  S  =  0,  a  is  given  by  (30)  with  t  instead  of 
t,  and  P  follows  from  (21).  At  the  free  boundary  on  the  right,  u  =  0, 

S  =  0,  P  =  0. 

The  stepsize  in  time,  AT,  is  determined  as  follows.  Of  the  two 
characteristics,  consider  the  one  which  has  the  larger  slope  with  respect  to 
the  t-axis  (for  example  the  one  labeled  I  in  Fig.  17);  such  a  slope  is  safely 
defined  by  |u|  +  a.  If  A,  B  and  C  are  the  three  points  used  in  the  compu¬ 
tation,  draw  a  line  issuing  from  C  and  symmetric  to  I.  Let  D  be  the  inter¬ 
section  of  I  and  its  symmetric.  Every  value  of  At  less  than  BD  satisfies 
the  Courant-Friedrichs-Lewy  rule  (Ref.  10).  Now,  BD  =  Ax/(u|  +  a)  and 
so  we  assume 

At  =  .7  — — —  (56) 

M  +  a 

The  coding  of  this  program  is  very  easy;  a  sample  code  is  shown  in 
Appendix  2  and  it  should  be  self-explanatory.  The  "output"  routine  is  not 
reported  here  and  is  left  to  the  reader  to  be  shaped  according  to  his 
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preferences.  The  "piston"  routine  is  not  shown  either;  it  must  provide 
the  necessary  data  to  define  the  piston  path,  b(t)  and  its  derivatives  b 

II 

and  b. 

2.  2  The  case  of  discontinuous  initial  acceleration. 

We  apply  now  the  program  to  the  two  problems  mentioned  in 
Section  2.  First,  we  assume  the  piston  path  as  defined  by  (41).  We  know 
that  u^  is  discontinuous  along  the  x  =  Jyt  line  and  that  it  becomes  infinite  at 

tA=7+T='4929  (for  Y  =  1.4)  (57) 

where  a  shock  should  originate.  We  are  interested  in  seeing  what  kind  of 

difficulties  and  inaccuracies  arise  because  of  the  discontinuities  in  the  first 

derivatives.  Such  discontinuities  are  not  as  strong  as  the  ones  across  a 

shock  but  still  they  make  the  flow  parameters  not  twice  differentiable  along 

the  first  I- characteristic.  So,  we  expect  a  progressive  worsening  of  the 

computed  flow  in  that  neighborhood,  until  a  pattern  similar  to  the  one 

described  in  Fig.  7  builds  up  as  representative  of  the  shock. 

$ 

Fig.  18  shows  some  of  the  results  (values  of  u)  .  The  computed 
values  are  shown  as  circles,  whereas  the  exact  distribution  is  represented 
by  a  continuous  line.  Point  A,  as  defined  by  (39),  is  also  shown.  We  see 
that,  long  before  the  characteristics  coalesce,  the  computed  values  deterio¬ 
rate  in  the  vicinity  of  the  line  where  the  derivatives  of  u  are  discontinuous. 

In  addition,  once  the  characteristics  coalesce  into  a  shock,  a  wiggle  starts 
building  up.  According  to  our  analysis  (Fig.  7),  we  should  fit  a  shock  two 

*  The  mesh  has  nodes  spaced  .025  apart.  There  are  no  other  nodes,  in  the 
range  of  interest,  except  those  shown  in  Fig.  3.18.  However,  not  all  the 
steps  are  shown.  The  computation,  with  a  total  of  50  intervals  across, 
takes  about  6  seconds  on  a  CDC  6600  computer  (including  compilation).  The 
computation  has  been  repeated  on  an  IBM  360/50  computer,  coupled  with  the 
Stromberg- Carlson  4020  plotter,  to  obtain  the  figure. 
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nodal  points  at  the  right  of  the  wiggle.  We  will  see  later  on  how  to  treat 
the  shock  as  u  discontinuity;  by  doing  so  (here  we  mention  the  result  with¬ 
out  presenting  any  figure),  the  shock  would  quickly  adjust  itself  to  the  right 
position,  and  the  values  on  the  high-pressure  side  would  improve  rapidly. 
However,  the  results  are  far  from  perfect  in  the  time  space  starting  at 
about  t  =  .4  and  lasting  for  a  time  of  the  order  of  1.  This  effect  is  due 
not  to  the  formation  of  the  shock  but  to  the  discontinuity  in  the  derivatives 
of  the  velocity.  We  will  return  on  the  subject  on  Section  3.3. 

2.  3  The  case  of  continuous  initial  acceleration. 

Now  the  same  program  is  applied  to  the  case  where  the  piston  path 
is  defined  by  (43).  The  derivatives  of  the  velocity  are  continuous  every¬ 
where  and  the  coalescence  of  characteristics  occurs  at  point  B,  defined  by 
(43).  With  y  =  1.4,  we  find  t  =  .389  and,  from  (38), 

Jd 

tQ  =  .844,  =  .8415  (58) 

Since  in  this  case  we  are  not  forcing  the  numerical  computation  to 
approximate  derivatives  where  they  do  not  exist,  we  anticipate  a  much 
better  behavior  of  the  computed  results  prior  to  point  B.  The  expectation 
is  confirmed  by  Fig.  19.  Again,  as  in  Fig.  18,  the  exact  distribution  of 
u(x)  is  represented  by  a  continuous  line,  and  the  computed  values  (at  mesh 
points  spaced  .2  apart)  by  circles.  The  circles  fall  exactly  on  the  lines, 
as  long  as  t  <  tg.  In  addition,  it  can  be  noted  that  the  numerical  technique 
is  capable  of  handling  a  very  sharp  transition  with  no  trouble. 

If  the  computation  is  continued  beyond  t  =  t_,  however,  the  effects 

1J 

of  neglecting  a  shock  begin  showing  up  and  a  wiggle  if  formed  where  antici¬ 
pated  (two  mesh  points  at  the  left  of  the  approximate  location  of  the  shock). 
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ARIOUS  TIMES  IN  THE  FLOW 
!  LAW  x„=t3 
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2.4.  An  expansion  problem  with  discontinuous  initial  acceleration. 

So  far,  we  have  found  experimental  confirmation  to  the  points  made 
above.  Discontinuities  in  the  flow  parameters,  such  as  shocks,  originate 
wiggles.  Discontinuities  in  their  derivathes  are  not  equally  catastrophic 
but  deteriorate  the  numerical  results  in  a  wider  and  wider  region  surround¬ 
ing  the  location  of  the  discontinuity.  In  both  cases,  the  discontinuities  should 
be  treated  as  boundary  points . 

Before  analyzing  how  this  can  be  done,  it  is  interesting  to  see  how 
the  discontinuity  in  the  derivatives  of  the  flow  variables  affects  an  expanding 
flow.  To  this  effect,  we  study  the  expanding  counterpart  of  the  case  pre¬ 
sented  in  Section  2.2.  Let  the  piston  now  move  to  the  left,  and  its  path  be 
defined  by 

x  =  -  t2  (59) 

Again,  (35)  and  (36)  show  that  u  and  u  are  discontinuous  on  the  first 

X  t 

I- characteristic.  For  example, 


u  = 


(y-i)t  +  v7 


(60) 


However,  ux  decreases  with  increasing  t,  and  we  may  anticipate  that  the 
effects  of  the  discontinuity  in  ux  will  become  less  and  less  relevant  with 
increasing  time;  and  since  they  are  not  strong  to  begin  with,  they  should  be 
irrelevant  throughout. 

The  anticipation  is  confirmed  by  Fig.  20,  drawn  with  the  same 
criteria  as  Figs.  18  and  19.  The  computed  values  (circles)  lie  very  close 
to  the  theoretical  line.  To  see  the  effects  of  neglecting  the  discontinuity 
one  has  to  enlarge  the  u- scale  by  one  order  of  magnitude. 


3.  Treatment  of  a  shock  on  a  boundary. 

In  this  section  we  examine  in  detail  a  technique  which  combines  the 
interior  point  computation  described  above  with  a  treatment  of  one  boundary 
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2.4.  An  expansion  problem  with  discontinuous  initial  acceleration. 

So  far,  we  have  found  experimental  confirmation  to  the  points  made 
above.  Discontinuities  in  the  flow  parameters,  such  as  shocks,  originate 
wiggles.  Discontinuities  in  their  derivatives  are  not  equally  catastrophic 
but  deteriorate  the  numerical  results  in  a  wider  and  wider  region  surround¬ 
ing  the  location  of  the  discontinuity.  In  both  cases,  the  discontinuities  should 
be  treated  as  boundary  points. 

Before  analyzing  how  this  can  be  done,  it  is  interesting  to  see  how 
the  discontinuity  in  the  derivatives  of  the  flow  variables  affects  an  expanding 
flow.  To  this  effect,  we  study  the  expanding  counterpart  of  the  case  pre¬ 
sented  in  Section  2.2.  Let  the  piston  now  move  to  the  left,  and  its  path  be 
defined  by 

x  =  -  t2  (59) 

Again,  (3  5)  and  (36)  show  that  u^  and  ut  are  discontinuous  on  the  first 
I- characteristic.  For  example, 

u  = - - - —  (6o) 

( y - 1 ) t  +  vV 

However,  u^  decreases  with  increasing  t,  and  we  may  anticipate  that  the 
effects  of  the  discontinuity  in  u^  will  become  less  and  less  relevant  with 
increasing  time;  and  since  they  are  not  strong  to  begin  with,  they  should  be 
irrelevant  throughout. 

The  anticipation  is  confirmed  by  Fig.  20,  drawn  ujfh  the  same 
criteria  as  Figs  18  and  19.  The  computed  values  (circles)  lie  very  close 
to  the  theoretical  line.  To  see  the  effects  of  neglecting  the  discontinuity 
one  has  to  enlarge  the  u-scale  by  one  order  of  magnitude. 

3.  Treatment  of  a  shock  on  a  boundary. 


In  this  section  we  examine  in  detail  a  technique  which  combines  the 
interior  point  computation  described  above  with  a  treatment  of  one  boundary 
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as  a  shock.  We  assume  that  the  location  of  the  shock  at  a  given  time  is 
known.  The  shock  is  moving  at  a  certain  speed,  W  =  c(t),  which  is  one 
of  the  unknowns  of  the  problem.  Consequently,  the  interior  point  computa¬ 
tion  must  also  be  modified  slightly,  and  we  begin  with  it. 


3.1.  Computation  of  interior  points. 

Referring  to  Fig.  16,  the  following  change  of  coordinates  is 
necessary: 

v  =  x  ~  W) 

c(t)  -  b(t) 

T  =  t 


(61) 


Since,  for  any  function  f, 

fx  =  ^Efx-  £,  =  £t+— <62> 

the  equations  of  motion  in  the  computational  space  are,  in  matrix  form. 


£T  +  A  £x  =  0 


with 


"p“ 

"C 

F 

0‘ 

f  = 

u 

,  A  = 

E 

C 

0 

-S- 

.0 

0 

c. 

(63) 


(64) 


i  r 

D  =  ,  C  =  D  j^u  +  (X-l)  b  -  Xc 


cj  ,  E  =  D7  ,  F  =  D’ 


Consequently, 


(65) 


Dx  =  °,  cx  =  D(ux+b- c*’  EX  =  D'7X*  FX~°*  AX  " 


CX  °  ° 

Ex  cx  0 
0  0  c 


and 


XJ 


(66) 


Dt=-D2(c),  CT=DTC/EH-D[uT+(X-l)b-XcJ,  Et=Dt7  +  d7t,  Ft=Dty. 


CT  FT  ° 
ET  CT  ° 


0 


0 
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The  code  then  proceeds  as  shown  in  Section  2.  1. 

3.2.  Shock  computation. 

To  compute  the  right  boundary,  where  the  shock  is,  we  proceed  as 
follows.  Since  the  flow  in  front  of  the  shock  is  known,  only  four  quantities 
must  be  evaluated  at  the  shock  point,  viz.  P,  S  (or  p  and  p ),  u,  and  W. 
The  Rankine-Hugoniot  conditions  can  be  used  here  in  the  form: 


=  wA/y 

(67) 

=  Fr<1  +  Ir-Mi> 

(68) 

_  (Y+l) n-(Y-l) 

“(Y+rr-(Y-rr0 

(69) 

V  =  -  W/o 

(70) 

u  =  v  +  W 

(71) 

where  p,  P»  v  and  u  are  values  on  the  high-pressure  side  of  the  shock, 
that  is,  at  the  point  we  wish  to  compute.  Here  we  have  five  equations  with 
six  unknowns  (W,  M^,  p,  p,  v.  u).  The  missing  equation  is  the  compati¬ 
bility  equation  on  the  I-characteristic  reaching  the  shock  point  from  a  point 
A  located  in  the  interior  of  the  region  being  computed.  In  the  physical  plane 


FIG.  21  COMPUTATION  OF  A  SHOCK  POINT 
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FIG.  22  NEWTON'S  RULE  TO  MINIMIZE  ERRORS 


the  situation  appears  as  in  Fig.  21 .  All  values  at  nodal  points  B,  C,  D,  ... 
are  known  at  time  t.  The  location  of  point  A  must  be  found  and  the  values  of 
the  physical  parameters  at  A  {all  denoted  by  *  in  what  follows)  must  be 
interpolated. 

To  determine  we  use  the  equation 

=  c  -  (u  +  a)  At  (?2) 


where  u  and  a  may  be  taken  as  the  arithmetic  averages  between  their 
values  at  Q  and  at  A.  The  compatibility  equation  (24-1)  can  be  written  in 

the  form: 


u 


ud 


“<P-  P*> 


(73) 


where,  again,  a  is  an  average  value. 

Equations  (72)  and  (73),  together  with  an  interpolation  scheme  to 
determine  u*.  P*.  a*,  must  be  added  to  the  Rankine-Hugoniot  equations 
above  and  the  whole  system  can  be  solved  by  iteration.  We  describe 


FIG.  23  EXAMPLES  OF  POSSIBLE  FAILURES  OF  NEWTON'S 
RULE 

Newton's  method  which  generally  works  in  a  variety  of  problems.  Two 
guesses  of  W  are  made,  Wj  and  W2>  Each  time,  equations  (67)  through 
(73)  are  used.  The  value  of  u  obtained  from  (73)  generally  is  different 
from  the  value  obtained  from  (71).  Call  E  their  difference  (Wj  when  Wj 
is  used,  E^  when  is  used).  An  extrapolated  value  of  W,  W,  which 
should  give  a  smaller  value  of  the  error  E,  is  obtained,  according  to 
Fig.  22,  as 


Then  the  procedure  is  applied  again,  using  instead  of  Wj,  W  instead 
of  and  instead  of  Ej,  and  reiterated  until  becomes  less  than 


'\  \| , ACTUAL  VALUE 


FIG.  24  INTERPOLATION  BEHIND  A  SHOCK 


a  prescribed  tolerance.  Such  a  technique  is  not  completely  safe.  For 
example,  there  are  possibilities  of  being  caught  in  a  loop  or  of  generating  a 
divergent  iteration  or  of  dividing  by  zero  even  when  a  solution  exists. 

Figure  23  shows  examples  of  these  possibilities  (the  exact  dependence  of 
E  on  W  being  shown  by  the  solid  line  and  the  solution  being  obviously  given 
by  point  S).  However,  such  cases  are  highly  improbable  because  they 
require  the  conditions  shown  in  the  figure  to  be  satisfied  to  an  order  of 
accuracy  equal  to  the  number  of  significant  figures  carried  by  the  computer. 
Therefore,  it  does  not  pay  to  take  special  provisions  for  them,  which  would 
make  the  code  cumbersome  and  increase  the  computational  time. 

More  interesting  because  of  its  physical  implications,  and  some¬ 
times  difficult,  is  the  problems  of  interpolating  values  at  A  (Fig.  21).  In 
most  cases  a  linear  interpolation  between  points  C  and  D  is  sufficient. 

This  amounts  to  approximate  the  value  of  any  function  f  at  A  by 


f*  =  f(C)  +  e 


f(D)  -  f(C) 


(75) 


x*  "  x(c) 


where 


There  are  cases,  however,  where  in  the  high  pressure  side  of  the 
shock  the  physical  values  are  distributed  as  shown  in  Fig.  24.  This  is  the 
case  at  the  very  beginning  of  the  shock,  in  both  problems  studied  above. 

If  point  A  is  very  clot  e  to  D  (and  this  happens  at  the  beginning  of  the  shock, 
when  W  is  practically  equal  to  u+a)  the  linearly  interpolated  values  at  A 
are  too  different  from  the  actual  values.  The  effect  of  the  pressure  surge 
behind  the  shock  is  not  picked  up  properly  by  the  computation  and  the  shock 
is  not  properly  fed.  Consequently,  the  strength  of  the  shock  does  not  grow 
sufficiently,  the  values  at  the  shock  do  not  grow  sufficiently  and,  in  few  steps, 
a  wiggle  starts  building  up  behind  the  shock.  The  effect  is  not  permanent 
since  the  right  values  are  transmitted  to  the  region  near  the  shock  from  the 
piston  and  the  flow  has  a  tendency  to  correct  itself.  See, for  example,  in 
Fig.  25,  the  u-distribution  in  the  problem  cf  Section  2.2  at  successive  times 
after  the  shock  formed,  if  the  technique  described  above  and  a  linear  inter¬ 
polation  at  point  A  are  used. 

A  much  better  numerical  description  of  the  physical  nature  of  the 
flow  in  the  vicinity  of  the  shock  is  obtained  if  the  values  at  A  are  considered 
as  the  average  of  the  values  extrapolated  linearly  from  points  B  and  C  and 
the  values  interpolated  linearly  between  points  C  and  D  .  This  means  that 
the  value  of  any  function  f  at  A  is  assumed  to  be 

f*  =  f(C)  +  e[f(D)  -  f(B)]/2  (77) 

The  results,  in  the  same  problem  shown  in  Fig.  25,  are  now  those  of 
Fig.  26.  The  shock,  as  one  can  see,  picks  up  strength  from  the  very 


#  This  is  indeed  a  very  easy  way  to  improve  the  interpolation.  One  could 
think  of  other  ways  of  doing  it,  for  example,  by  interpolating  on  second 
order  fits  of  x(u),  x(P)  and  x(S)  but  such  interpolations  would  reauire  a 
great  deal  of  additional  computations  which  seems  to  be  unjustified,  in  the 
light  of  the  results  obtained  by  using  (77). 


FIG.  25  COMPUTED  DISTRIBUTIONS  OF  u(x)  AT  VARIOUS  TIMES  IN  THE 
PRODUCED  BY  A  PISTON  ACCELERATED  BY  THE  LAW  xD=t2 
CONSTANT  SPEED  IS  REACHED  (FIRST)  w 
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FIG.  26  COMPUTED  DISTRIBUTIONS  OF  u(x)  AT  VARIOUS  TIMES  IN  THE  FLOW 
PRODUCED  BY  A  PISTON  ACCELERATED  BY  THE  LAW  xn  =  t2  UNTIL 
A  CONSTANT  SPEED  IS  REACHED  (SECOND)  P 


FIG.  32  COMPUTED  DISTRIBUTIONS  OF  u(x)  AT  VARIOUS  TIMES 
FOR  THE  FLOW  PRODUCED  BY  A  PISTON  MOVING  WITH 
THE  LAW  Xr\ s  t3 
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beginning  and  no  wiggles  appear. 

Obviously,  if  the  distributio  of  values  in  the  vicinity  of  the  shock  is 
not  as  critical  as  in  Fig.  24,  (77)  is  numerically  equivalent  to  (75). 

Once  the  values  at  the  shock  are  determined,  including  W,  the 
shock  location  at  t  +  At  is  obtained  by  simply  writing  that 

c(t  +  At)  =  c(t)  +  W(t)  At  (78) 

3 . 3  The  infinitely  weak  shock. 

Being  able  to  compute  a  shock  point  on  the  boundary,  we  have  the 
case  of  the  accelerating  piston  with  a  finite  initial  acceleration  under  control. 
The  proper  way  of  computing  it  consists  of  assuming  that  a  shock  exists  from 
the  beginning,  issuing  from  the  origin  of  the  (x,t)  plane  with  a  speed  equal 
to  Vy  .  In  so  doing,  the  discontinuity  in  u^  does  not  affect  the  computation 
since  it  takes  place  rt  a  boundary  point,  which  is  assumed  to  be  a  locus  of  dis¬ 
continuities  and  is  treated  as  sucn.  On  the  other  hand,  an  infinitely  weak 
shock  is  a  characteristic  and,  since  the  assumed  shock  does  not  exist  before 
t  =  t^,  during  this  time  the  boundary  is  the  first  I- characteristic  itself. 

Point  A  (Fig.  10)  is  then  reached  as  a  boundary  point  and  the  computation  can 
be  continued  beyond  it  indefinitely  since  all  the  computational  apparatus  to 
take  care  of  a  growing  shock  is  available.  The  results  (Fig.  26)  are  remark¬ 
ably  accurate  and  the  transition  from  the  simple  wave  to  the  shocked,  non- 
isentropic  flow  is  extremely  smooth. 


In  Fig,  25  and  26  the  piston  path  is  defined  by 

b_  t2  (0<  .t  <  .  7395  ) 

1.  479t-.  5469  (t  >  .  7395) 

By  doing  so,  a  steady  state  should  be  reached  asymptotically,  where 
the  shock  Mach  number  is  2  and  the  flow  is  uniform  between  the  piston 
and  the  shock.  It  can  be  seen  from  both  figures  that  such  a  steady 
state  is  reached  very  rapidly.  The  computed  shock  Mach  number  at 
the  time  corresponding  to  the  last  plot  is  2,01, 
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One  word  is  in  order  about  how  to  start  the  computation.  Since  the 
physical  region  to  be  computed  vanishes  at  the  origin,  we  are  forced  to  start 

at  some  very ‘small  but  finite,  positive  time.  In  general,  when  a  similar 
situation  occurs  in  a  more  complicated  problem,  one  may  assume  that  the 
initial  distribution  of  physical  parameters  corresponds  to  a  simple  wave. 
Appendix  3  contains  the  code  for  the  computation  from  which  Fig.  26 

was  obtained.  The  subroutines  "piston"  and  "output"  are  omitted.  By  com¬ 
paring  Appendix  3  with  Appendix  2  one  can  see  that  the  subroutine  which 

computes  the  interior  points  is  not  made  more  difficult  by  the  stretching  of 
the  physical  space  into  a  computational  space. 

A  further  proof  of  the  possibility  of  considering  a  character¬ 
istic  as  an  infinitely  weak  shock  is  found  if  the  code  of  Appendix  3  is 
applied  to  the  expansion  problem  of  Section  2.  4.  The  results  are  omitted 
here.  Starting  with  16  mesh  intervals,  values  correct  to  within  5  signifi¬ 
cant  figures  at  least  are  obtained  as  long  as  the  mesh  size  does  not  become 
gigantic. 

3.4  Treatment  of  an  imbedded  shock. 

In  the  problem  of  an  accelerating  piston  with  zero  initial  accelera¬ 
tion  the  shock  is  imbedded  in  the  flow  field.  In  this  section  we  describe 
briefly  how  this  situation  is  handled. 

Since  the  inception  of  the  shock,  that  is,  since  t  =  t_,,  the  flow  field 

D 

is  split  into  two  regions,  bounded  by  the  shock.  There  are,  thus,  two  points 
representing  the  shock,  one  being  the  right  boundary  of  the  high  pressure 
region  and  the  other  being  the  left  boundary  of  the  low  pressure  region 
(Fig.  27).  The  computation  of  the  former  is  the  same  as  outlined  in  Section 
3.2,  once  the  low  pressure  point  has  been  evaluated.  The  latter  is  reached 
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FIG. 27  TWO  COMPUTATIONAL  REGIONS 
SEPARATED  BY  A  SHOCK 


FIG.  28  COMPUTATION  OF  THE 
LOW-PRESSURE  SIDE 
OF  A  SHOCK 


FIG.  29  COMPUTATIONAL  REGIONS 
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by  two  characteristics  and  a  particle  path;  therefore,  two  compatibility  con¬ 
ditions  and  the  condition  of  conservation  of  entropy  can  be  written  to  deter¬ 
mine  the  flow  parameters  (Fig.  28).  In  the  physical  plane,  the  locations  of 
points  Aj,  and  B  are  given  by 

X*1  =  XQ  "  ^U+a*  At 

x#2  =  XQ  ”  (u-a)  At  (79) 

x  =  x  -  u  At 
V  B  Q 

After  the  pertinent  interpolations  (in  the  spirit  of  (77)  )  have  been  made, 

equations  (24)  are  applied,  in  the  approximate  form: 

p  _  a*l  P*1  +  a*2  P#2  +  y^u*2  ‘  u*l^ 
a*l  +  a*2 

u  =  u*i  - -r  (p- p*i»  <80> 

S  SB 

The  interior  points  in  each  region  are  treated  as  explained  in 
Section  3.1,  with  obvious  modifications  (for  example,  what  is  c(t)  for  region 

1  is  b(t)  for  region  2). 

The  right  arbitrary  boundary  of  region  2  can  be  chosen  to  be  the 
x  =  VY  t  line,  as  in  the  problem  with  finite  initial  acceleration.  Values  along 
that  line  can  be  safely  forced  to  be  u  =  0,  P  =  0,  S  =  0.  The  line  itself  is 
prescribed,  so  that  in  region  2,  c  =  t  ,  c  =  >Jy,  c  =  0. 

There  is  a  time  t^  (Fig.  29)  when  the  shock  between  regions  1  and 

2  crosses  the  first  I-characteristic.  In  its  vicinity,  if  no  special  provisions 
are  taken,  the  time  step  would  decrease  indefinitely  since  it  is  proportional 
to  tile  physical  mesh  size  and  this  in  turn  decreases  indefinitely,  when  the 
region  2  is  squeezed  between  the  shock  and  the  first  characteristic  (  the 
number  of  mesh  points  remaining  constant).  Therefore,  once  the  width  of 
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FIG.  30  BEHAVIOR  OF  ulx)  IN  THE 
OF  AN  INCIPIENT  SHOCK 
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region  2  becomes  less  thana  prescribed  value,  the  computation  of  the  region 
is  dropped  and  only  region  1  is  computed,  letting  the  shock  divide  it  from  the 
gas  at  rest. 

3.5.  Prediction  of  the  formation  of  a  shock. 

At  t  =  t^,  one  of  the  two  regions  disappears.  The  inverse  situation 

c 

arises  at  t  =  tg,  and  it  is  a  much  more  difficult  one.  Prior  to  t^,  the 
computation  proceeds  as  described  in  Section  2.  3.  We  need  a  criterion  to 
detect  the  origin  of  a  shock,  and  the  criterion  must  furnish  the  location  of 

point  B,  in  space  and  time,  to  within  one  mesh  interval  in  space  and  one 
or  two  steps  in  time,  lest  wiggles  appear.  This  is  indeed  a  very  difficult 
problem  in  that  a  point  must  be  found  where  u^  becomes  infinite,  when  we 
only  dispose  of  a  discreet  number  of  values  of  u.  Obviously,  no  finite 
difference  approximation  to  a  derivative  can  ever  become  infinite  and  any 
identification  of  infinity  with  "a  very  large  value"  is  arbitrary.  A  study  of 
the  flow  in  the  vicinity  of  point  B  (Fig.  19)  should  guide  us  to  find  a  better 
criterion. 

Let  us  consider  the  mesh  interval  where  the  slope  of  u(x)  is  the 
highest  and  four  nodes  (in  the  computational  space),  two  at  the  left  and  two 
at  the  right  of  it  (Fig.  30).  Since  we  are  interested  in  a  situation  where 
ux  becomes  infinite,  it  is  better  to  consider  x  as  a  function  of  u  rather 
than  u  as  a  function  of  x. 

The  simplest  polynomial  having  a  chance  of  giving  a  fair  approxima¬ 
tion  to  x(u)  in  the  interval  shown  in  the  figure  is  a  third  order  polynomial: 

x  =  Au^  +  Bu^  +  Cu+D  (81) 

We  are  looking  for  the  case  where  the  condition 

^=3Au2+2Bu+C  =  0  (82) 

is  satisfied  by  one  real  value  of  u.  This  requires 
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B2  -  3  A  C  =  0  (83) 

Suppose  now  that  A,  B,  and  C  are  determined  to  fit  the  values  of 

x(u)  at  the  four  points  of  Fig.  30.  We  can  use  (83)  as  a  test.  In  the  simple 
2 

wave,  B  -  3 AC  should  generally  be  negative  since  no  real  value  of  u 

2 

satisfies  (82).  When  t_,  is  reached,  B  -  3AC  should  vanish  and  then 

-D 

become  positive,  and  this  yields  a  possible  criterion  for  the  detection  of  tg. 
If  Xj  is  the  value  of  x  at  point  1,  the  values  of  x  at  points  2,  3,  and  4  are 
Xj  +  Ax,  Xj  +  2  Ax,  Xj  +  3  Ax.  By  writing  (81)  for  these  four  points,  and 
calling  6  the  determinant  of  the  coefficients  of  the  right  hand  sides,  we 
obtain 


A 


_  Ax 

“  T 


0  Uj  Uj 


1  u2  u2 


2  u3  u3 


3  u4  u4 


1 

1 

1 

1 


Ax 

T 


U 


U, 


3 

2 


+  u; 


+  u. 


u3-2u2  +  ui 

u4-2u3  +  u2 


and,  similarly, 

B  =  ■y-  f ( 3 , 1 )  ,  C  =  -  ^f(3,2) 


Therefore  the  criterion  suggested  by  (83)  can  be  written  as 

[f(3, 1 ) J2  -  3  f ( 3 , 2 )  £(2,1)  <  0 


(84) 


Fig.  31  shows  the  trend  of  the  left  hand  side  of  (84)  in  the  calculation  of 
Fig.  19. 

Finally,  Fig.  32  shows  plots  of  u(x)  at  different  times  for  a  calcula¬ 
tion  of  the  problem  of  section  2. 3  carried  along  indefinitely  beyond  tn.  The 

D 

criterion  above  to  determine  tg  is  applied,  then  the  computational  region  is 
split  into  two  parts;  at  t^  the  second  region  is  dropped  again,  as  described 
above.  The  coding  of  this  program  is  not  shown  since  it  can  be  easily 
obtained  as  an  extension  of  Appendix  3.  All  simple  arrays  should  in  this  case 
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be  double  arrays,  to  permit  values  of  the  physical  parameters  in 
different  regions  to  be  distinguished  from  each  other;  similarly, 
quantities  such  as  b,  c,  b,  tsc  etc.  ,  should  be  stored  in  arrays,  one 
value  for  each  region.  More  details  on  the  treatment  of  a  problem 
with  multiple  shocks,  including  the  problem  above  as  a  particular 
case,  will  be  found  in  another  paper. 
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APPENDIX  I 


DIMFNSION  R  <  50 ) .0(50) ,E(50) .QR(50) #QM(50) ,HR(50) #HM (50 ) ♦ HE (50 ) *RN 
1  (50)  ,EN(50)  *U(50)  »H(e;n) 

REAL  M  ( 50 ) » MN ( 50 ) 

DATA  SIG.P?.R2»E1.Q1,0E/. 1,4. 5. 2. 66666666666667. 5. 3.6. 6.. 4/ 

100  FORMAT ( 1H0/ ( 10E12. 4) ) 

101  FORMAT ( 1H1/I10) 

SIG?sSlG*#?  S  GF=SQRT (1.4)  I  U1=-2.*GF 

U2=IJ1/R2  5  E?=2.5*P2  +  1.05  S  H1*-12.6*GF 

DOl  m=1*20  S  F(N)=E2  S  U(N)«U2 

1  R  ( N) =R2 

DO  2  N=21.40  $  E(N)=E1  S  U(N)*U1 

2  R (N) si  • 

DO  3  N=1.40  S  M(N)*U1  $  Q(N)=Q1 

3  H  ( N ) sH 1 

U  ( 2 1 ) =Ul*.000l 
M ( 2 1 ) =U ( 2 1 ) 

0(21 ) ».4*E (21 ) ♦.8*U(21 ) **2 

H(2i ) =U (21 ) * ( 1  .4*E ( 2 1 )-.2*U(21)**2) 

KsO 

8  K  =  K  ♦  1 

IF (FLOAT (K/l 0 ) .NE .FLOAT (K ) /l 0 . ) GO  TO  7 

WRITF(6.101)K 

WRlTF (6.100) (U  (N) *N*1 .40 ) 

*/RI TF  (6.100)  (M(N)  »N=1.40) 

WRITF (6.100) (Q  (N) »N=1 *40 ) 

7  IF (K.GT.200) CALL  EXTT 
DO  4  N= 1.40 
QR(M)=-.8*U(M>**2 
QM(N)  =1  .b*IJ  (N) 

NR(M)=-1.4»U(N)*E(N)/R(N)^.4*U(N)**3 
HM(N)=1.4*F (N)/R(N)-.A*U(N)#*2 

4  HE(M)=1.4*  (M) 

UO  <=.  N=2.39 

DMPsm(N*1)-M(N)  $  DQP=Q(N^1 ) -Q(N)  5  DHP*H (N+l ) -H (N) 
DMMsm(N)-M(N-1)  S  DQM=Q(N)-Q(N-1)  S  DHMsH(N)-H(N-1) 

RN(N) =R (N) - (OMP.DMM) «SI6*2.» (DQP-DQM)*SIG2 

MN(N)=M(N)-(DQP*DQM)*SIG* ( ( OR ( N. 1 ) ♦OR (N) ) *DHP* ( OR ( N) ♦QR (N-l ) ) *DmM* 
1  (UM(nO) ♦Om(N) )*DQPMQM(N) ♦QM(N-i) ) *DQM^2,*QE* (DHP^DHM) )*SlG2 
EN (N) SE ( N) -  (DHP  +  DHM)  *SlG  +  ( HR  ( N«  1 )  ♦  HR  ( N >  >  *DMP+  (HR  (N)  ♦HR(N-l)  )#DMM* 
1 (HM<mO) ♦HM(N) )*DQP*(HH(N)4HM(N-1) ) *DQM  + (HE (N+ 1 ) *HE (N) )*DHP.(HE(N) 
2*HE ( M— 1 ) )*0HM)*SIG2 
IF ( APS (RN (N) -R (N) ) .LT.l.E-R)  RN(N)=R(N) 

I F ( APS ( MN ( M ) -M ( N ) ) .Lt.I.E-P)  MN(N)sM(N) 

IF(ARS(EN(N)-E(N) ).LT.1.E-o)  EN(N)=E(N) 

5  CONTINUE 

DO  6  N=2  * 39  $  R(N)=RN(N)  $  M  (N) *MN (N) 

E (N) sEN (N)  S  U (N) =M (N) /R (N)  S  Q (N )  =  .  4«E ( N) ♦ . 8*U (N) *M (N) 

6  H(N)«U(N)*(1.4*E(N)-.?*U(N)*M(N) ) 

GO  TO  8 

END 


I 
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APPENDIX  II 


COMMON  K  tNA  *NC » GAMMA  txOtTENDtGAtQDtGFtTlMFtDT  t DT2  t DX  t DDX t DDX2 1 BSEC 
1 «  X ( 1 O  0 ) t  XX  ( 1 0 0 ) *  P ( 1 0  0 )  *  U ( I o  0 )  *  A ( 1 00 ) *  T ( 1 0  0 ) *S(100) » SN ( 1 00 ) » UN ( 1 0  0 ) 
2*PN ( 1 00 ) 

DATA/X/1000*0./ 


100  FORMAT (HE10. 4) 

101  F0RmaT(1H1,4F15.4) 

REA!)  (5«  100)  GAMMA tX0t AN t  TEND 
WRlTF ( 6 « 1  0 1 ) GAMMA  «XO»aN»TENO 
GA=OaMMa/ (OAmMA-1 • ) 

GD=.*i*  (GAMMA-1 . ) 

GFsSORT (GAMMA) 

OX--XO/AN 

NA  =  AM 

UDX=.b/DX 

DDX?rl ,/UX/DX 

NC  s  M  A  ♦  1 

DO  ?  N=1,NC 

X (N) =DX*FLOAT (N-l ) 

XX  <  N ) =  X (N) 

A (N) =GF 

2  T  ( N )  s  1  • 

K=1 

TIMFrO. 

3  KsK ♦  1 

0T*100. 

DO  7  N=1»NC 
0T1=.7*UX/(U(N)4A(N) ) 
IF(OTl.GT.DT)GO  TO  7  ■ 

DT=OT 1 

7  CONTINUE 
TIMFrTlME^nT 
UT2rpT**2/2 • 

CALL  PISTON 
Al=GF*GDttUN(l) 

PN( 1 ) =GA*Al0G(A1**2/GAMMA) 
CALL  POINTS 
DO  R  N=1»NC 
XX(M)«XX(1)*X(N) 

U (N) -UN (N) 

S(N) rSN(N) 

P(N) =PN(N) 

T (N) sEXP (P (N) /GAtS (N) /GAMMA) 
A(N) =GF*SUrT(T(N) ) 

8  CONTINUE 
CALL  OUTPUT 

IF(TIME.LE.TEND)G0  TO  3 

CALL  EXIT 

END 


subroutine  POINTS 
DO  1  N=2  tNA 

UX=(t)(Nol)  -U(N-1  )  )*DDX 
px=(p(n+d-P(n-1)  )*Dnx 
SX=(S(N^1)-S(N-1) ) *DDx 
TX=T(N) * (PX/GA*SX/GAmmA) 

UXXs<U(N+l)  ♦IJ(N-1)-2,*U(N)  )*DDX2 
PXX={P(N+1)*P(N-1)-2,*P(N) ' *DDX2 
SXXs(S(N*i) ♦S(N-1)-2,*S(N) ) *DDX2 
AC=U ( N ) -U ( 1 > 

PT=-(AC*PX-GAMMA*UX) 

UT=- (AC*UX*T (N) *PX) 

ST=-aC*SX 

TT=T (N) # (PT/GA* ST/GAMMA) 

ACXsiiX 

ACT*UT-BSEC 

PTXa-  (ACX*PX*AC*PXX*GAMMA*UXX) 

UTXs-  (ACX*I)X*AC*UXX*TX*PX*T  (N)  *PXX) 
STX=- (ACX*SX*AC*SXX) 

PTT=- (ACT*PX*AC*PTX*GAMMA*UTX) 
UTTs-(ACT*UX*AC*UTX*TT*PX*T (N) *PTX) 
STTs- (ACT*SX*AC#STX) 

PN(N) =P (N) ♦PT#OT  tPTT  *OT2 
UN (N)3 *  S * 7 8U (N) ♦UT*OT*UTT#qT2 
SN(M)=S(N) ♦ST*UT*STT*DT2 
1  CONTINUE 
RETURN 
END 
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APPENDIX  III 


COMMON  K»NA*NC* GAMMA, x0» TEND *GA,SB»GC»G0?GE *GF* TIME »DT *DT2»D*DX, 

1 OOX, D0X2 td,C, BOOT. COoTfBSECtCSEC»XA.X( 100) ? XX ( 1 00 ) »P ( 1 00 ) »U ( 1 00 ) » 
2A (100) »T(100) .S(100) »PN(100) *UN(100) »SN(100) 


c 

100  format (8E10. 4) 

101  FOHmaT ( 1H1 ,4E15.4) 

C 

REAR (5* 100) GAMMA, XO* AN* TEND 
WRITE (6* 101 ) GAMMA *X0* AN* TEND 
GA=r,AMMA/  (GAmMA-1  .  ) 

GH  = 1 ./ (GAMMA-l . ) 

GC= ( GAMMA* 1 • ) *GB 
GD=,S* (GAMmA-1 . ) 

GE=,«;*  (GAMMA*1.  ) 

GF=SORT (GAMMA ) 

N  A  =  A  N 
NC  =  MA*  1 
TIMFr.Ol 
CALL  RISTON 
H  =  X  X  (  1  ) 

U(l)aUNd) 

XX (NC ) =GF*TImE 
DXs  1  ./AM 
0=XX (NC) -XX ( 1 ) 

DO  X= . S/l)X 
DDXPsl ./DX/DX 
DO  ?  N= 1 » Mr 
X (M) =DX*FLOAT (N-l ) 

XX (M) =X ( N) *D*XX ( 1 ) 

U ( N ) rU(l)*(l,-X(N) ) 

P(N) =2.*GA*AL0G(1.^U(M) *GD/GF) 
S(N) -0. 

A (N) =GF ♦GD*U(N) 

2  T (N) =A (N) «*2/GAMMA 
C=XX (NC) 

CDOT  =  GF 
CSEC=0. 

K  =  1 

3  KsK ♦ 1 
DT=1 no. 

DO  7  N=1«NC 

0T1=,7#UX/{U{N)*A(N) ) *D 
IF  (DTI .GT.DT)GO  TO  7 
UT=Ot1 
7  CONTINUE 

timf=time*ot 

UT2=0T«*2/2. 

CALL  PISTON 
Al=GF*GO»UN{ 1 ) 

PN< 1 ) =GA*ALOG (A1**2/GAMMA) 

SN(n  =0. 


CALL  points 
CALL  SHOCK 
BsXX  ( 1 ) 

BDOT=UN(l) 

D=C-H 

DO  A  N=1*NC 
XX (N) =B*X (N) *D 
U(N) =UN(N) 

P(N)=PN(N) 

S(N)sSN(n> 

T (N) sEXP (P (N) /GA*S (N) /GAMMA) 
A(N)rGF*SQRT(T(N) ) 

8  CONTINUE 
CALL  OUTPUT 

IF (TIME.LE.TEND)GO  TO  3 

CALL  EXIT 

END 
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APPENDIX  III  (con1 1.  ) 


SUBROUTINE  points 
DO  1  N=2*NA 
AAsl ,/D 

UX=(U<N*1)-U(N“1) )*DDX 
PX=(P(N*1)-P{N-1) )*DDX 
SX=(S(N*1)-S(N-1))*DDX 
TX=T (N) * (PX/GA*SX/GAmmA) 

UXXs  (U(N*1)  ♦IJ(N-1)-2.*U(N)  ) *DDX2 
PXXr (P (N*l ) ♦P(N-1)-2.#P(N) ) *DDX2 
SXX*(S(N*1)*S(N-1)-2.#S(N) )*DDX2 
AC=AA*( <X(N)-1.)*BD0T-X(N)*CD0T*U(N) 
AD=GAMMA*AA 
AE*AA*T (N) 

RTs- ( AC*PX*AD*UX) 

UT=- ( AC"UX*AE*PX) 

ST  =  -AOSX 

TT  =  T (N) * (PT/GA+ST/GAMMA) 

ACXrAA*(dDoT-CDOT*UX> 

AEXsAA*TX 

PTXs- (ACX*PX*AC*PXX*AD*UXX) 

UTX=- (ACX*UX*AC*UXX*AFX*PX*AE*PXX) 

STX=- (ACX*SX*AC*SXX) 

AAT=-AA**2* (CCOT-BDOT) 

ACT=AAT»AC/AA*AA*{ (X(N)-l. ) *BSEC-X (N)*CSEC*UT) 

ADTrOAMMA*AAT 

AET* A AT*T ( N ) *AA*TT 

PTT  =  -  (ACT*pX*AC*PTX*AOT*UX<-AD*UTX) 

UTTs- (ACT*UX*AC*UTX>aFT*PX*AE*PTX) 

STT=- (ACT*SX*AC*STX) 

PN(N) =P (N) *PT*UT*PTT*OT2 
UN(N)=U(NUUT*DT*UTT#DT2 
SN<N)=S(N)  ♦ST*DT*STT#OT2 
I  CONTINUE 
RETURN 
END 
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APPENDIX  III  (con't.  ) 


SUBROUTINE  SHOCK 
DIMENSION  W<2) t ERR (2) 

102  FORMAT (8H0FAILURE) 

C=CDOT*DT*C 

ME  =  1 

KIPal 

W  ( 1 ) aCDOT 

W  ( 2 ) =1.000l«CDOT 

UN(NC) =U(NC) 

AN=A (NC) 

UST  AD  =  U { NC ) 

ASTAP=A (NC) 

60  EM?*w (ME) **2/GAMMA 
PS=(GAMMA*FM?-GO)/GE 
RS=  («C*PS*1  .  )  /  (GC+PS) 

XXSTAR=C- (UN(NC) ♦USTaR*AN*ASTAR)/2.*DT 
XSTAR= (XXSTAR-B) /d 
EPSs (XSTAR-X (NA) ) /OX*. 5 
ASTAP=A(NA) ♦EPS*(A(NC)-A(NA-1) ) 

USTAR=U(NA) *EPS*(U(NC)-U(NA-1) ) 

PSTAP=P(NA) ♦FPS*(P(NC)-P(NA-i) ) 

UREl  =-W  (ME  )  /PS 
UN  (NO  )  =UREl_*w  (ME  ) 

PN(NC) =PSTAR-GAMMA/ASTAR* (UN(NC)-USTAR) 
AN=SoRT (GAMMA*PS/RS) 

ERR(mE)PN(NC)-ALOG(PS) 

IF ( MF • EQ. 2 ) GO  TO  62 
ME  =  ? 

GO  TO  60 

62  IF(ARS(ERR(ME) ) .LE.1.E-5)G0  TO  61 

WA=w ( 1 ) -ERR ( 1 ) * (W (2) -W(l ) ) / (ERR (2) -ERR ( 1 ) ) 
ERR ( i  ) =ERR ( 2 ) 

W ( 1 ) sW (2) 

W (2) sWA 
KIPsKlP*l 

IF(KIP.LE.20)GO  TO  60 
WRITF ( 6 « 10?) 

CALL  EXIT 

61  SN(NC)=PN(NC)-ALOG(RS)*GAMMA 
CSEC= (W(ME)-CDOT)/DT 
COOTrW (ME) 

RETURN 

END 
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